#include "SparseMatrix.h"


/**
 * the number of nnz for each row is stored in the first element
 *  col_idx[n_max_per_row*i], the first nnz is in col_idx[n_max_per_row*i + 1]
 *  we do not check the validaty of row and col, 
 */
/*
void sparse_add_real(int m_per_row, int *col_idx, double *val_nnz, 
		     int row, int col, double val)
{
  int pos;
  int begin = m_per_row*row;
  int nnz_row = col_idx[begin];
  for(pos = begin + 1; pos < begin + nnz_row + 1; pos++){
    if(col_idx[pos] == col);
  }
  col_idx[pos] = col;  // if exist, keep unchange
  val_nnz[pos] += val;
  col_idx[begin] = (pos > begin + nnz_row)?pos:nnz_row;
}
*/
/**
 * treat the sparsity pattern from matrix J0 and S0 the row compress
 * storage format is gotten : I,J,S we assume the memory is reserved by
 * calling functions.  for ease of matrix iteration method, the diagnal
 * entries are put at the beginning of each row.
 */
 /*void sparse_compress_real(int m_per_row, int *col_jdx, double *val_nnz, 
			  int *row_ptr, int *col_idx, double *val, int dim)
{
  int i,j,k,pos,idx;
  // find out start position of each row
  row_ptr[0] = 0;
  for(i = 0; i < dim; i++){
    row_ptr[i+1] = row_ptr[i] + col_jdx[m_per_row*i];
  }
  // then append the nnz
  for(i = 0; i < dim; i++){
    idx = row_ptr[i] + 1;  // researve the first place for diagnal element
    for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++){

      // the j'th element in col_jdx and val_nnz are corresponsing to
      // the k'th element in col_idx and val
      j = m_per_row*i + 1 + (pos - row_ptr[i]);
      k = idx;
      if(col_jdx[j] == i){
	k = row_ptr[i]      // remember that pos=0 is for the diagnal element
      }
      else{
	idx ++;
      }

      col_idx[k] = col_jdx[j];
      val[k] = val_nnz[j];
    } 
  }
}
 */

void sparse_pattern_Aphi(int **row_ptr, int **col_idx, int *nnz,
			int *Element,int *ElementEdge,
			int N_node, int N_edge,int N_element)
{
  int Msize = N_edge + N_node;

  int i = 0, j = 0, elem = 0;
  int pos = 0, row_start = 0;
  int exist = 0;

  // 1. malloc local variable for the calculation of the sparsity 
  int nz_max_per_row = 250;

  int *row_nnz = (int *)malloc(Msize*sizeof(int));
  int *col_ldx = (int *)malloc(Msize*nz_max_per_row*sizeof(int));
  int elem_dof[10];

  for(elem = 0; elem < N_element; elem++){

      // archive the degree of freedoms for current element
      for(i = 0; i < 6; i++) elem_dof[i]   = ElementEdge[6*elem + i];
      for(j = 0; j < 4; j++) elem_dof[6+j] = N_edge + Element[4*elem + j];

      // add each entries to global sparsity pattern
      for(i = 0; i < 10; i++){
	row_start = nz_max_per_row*elem_dof[i];
	for(j = 0; j < 10; j++){
	  
	  exist = 0; 
	  for(pos = row_start; pos < row_start + row_nnz[elem_dof[i]]; pos++){
	    if(col_ldx[pos] == elem_dof[j]){
	      exist = 1;
	      break;
	    }
	  }

	  if(!exist){  // do not exist yet, then add it!
	    col_ldx[pos] = elem_dof[j];
	    row_nnz[elem_dof[i]] += 1;
	  }

	} // finish 10 cols
      }  // finish 10 rows

  } // run over all element
  
  // 2. compress the sparsity pattern
  *nnz = 0;
  for(i = 0; i < Msize; i++){
    if(row_nnz[i] >= 250){
      printf("Too small max_nnz_per_row!\n");
    }
    (*nnz) += row_nnz[i];
  }

  (*col_idx) = (int *)malloc((*nnz)*sizeof(int));
  for(i = 0; i < (*nnz); i++) (*col_idx)[i] = -1;
  (*row_ptr) = (int *)malloc((Msize+1)*sizeof(int));
  for(i = 0; i < Msize+1; i++) (*row_ptr)[i] = -1;

  // compress the sparse pattern
  (*row_ptr)[0] = 0;
  for(i = 0; i < Msize; i++){
      for(j = 0; j < row_nnz[i]; j++)
	(*col_idx)[(*row_ptr)[i] + j] = col_ldx[nz_max_per_row*i + j]; 
      (*row_ptr)[i+1] = (*row_ptr)[i] + row_nnz[i]; 
    }

  // 3. free local variable
  free(row_nnz);
  free(col_ldx);
}


   /**
    * calculate sparse pattern for nodal finite element space
    * 
    */

void sparse_pattern_nodal(int **row_ptr, int **col_idx, int *nnz,
			  int *Element, int N_node, int N_element)
{
  int Msize = N_node;

  int i = 0, j = 0, elem = 0;
  int pos = 0, row_start = 0;
  int exist = 0;

  // 1. malloc local variable for the calculation of the sparsity 
  int nz_max_per_row = 250;

  int *col_ldx = (int *)malloc(Msize*nz_max_per_row*sizeof(int));
  int *row_nnz = (int *)malloc(Msize               *sizeof(int));
  for(i = 0; i < Msize; i++) row_nnz[i] = 0;

  int dof_i, dof_j;
  for(elem = 0; elem < N_element; elem++){

    for(i = 0; i < 4; i++){
      dof_i = Element[4*elem + i];
      row_start = nz_max_per_row*dof_i;
      for(j = 0; j < 4; j++){
	dof_j = Element[4*elem + j];

	exist = 0;   // search the entry pos
	for(pos = row_start; pos < row_start + row_nnz[dof_i]; pos++){
	  if(col_ldx[pos] == dof_j){
	    exist = 1;
	    break;
	  }
	}

	if(!exist){  // do not exist yet, then add it!
	  col_ldx[pos] = dof_j;
	  row_nnz[dof_i] += 1;
	}

      } // finish 4 cols
    } // finish 4 rows

  }

  // 2. compress the sparsity pattern
  *nnz = 0;
  for(i = 0; i < Msize; i++){
    if(row_nnz[i] >= 250){
      printf("Too small max_nnz_per_row!\n");
    }
    (*nnz) += row_nnz[i];
  }

  // initial the return variable
  (*col_idx) = (int *)malloc((*nnz)*sizeof(int));
  for(i = 0; i < (*nnz); i++) (*col_idx)[i] = -1;
  (*row_ptr) = (int *)malloc((Msize+1)*sizeof(int));
  for(i = 0; i < Msize+1; i++) (*row_ptr)[i] = -1;

  // compress the sparse pattern
  (*row_ptr)[0] = 0;
  for(i = 0; i < Msize; i++){
    for(j = 0; j < row_nnz[i]; j++)
      (*col_idx)[(*row_ptr)[i] + j] = col_ldx[nz_max_per_row*i + j]; 
    (*row_ptr)[i+1] = (*row_ptr)[i] + row_nnz[i]; 
  }

  // 3. free local variable
  free(row_nnz);
  free(col_ldx);
}


/**
 * Sparse pattern calculation for Edge element(P1)
 * 
 */
void sparse_pattern_edge(int **row_ptr, int **col_idx, int *nnz,
			 int *ElementEdge, int N_edge, int N_element)
{
  int Msize = N_edge;
  int i = 0, j = 0, elem = 0;
  int pos = 0, row_start = 0;
  int exist = 0;

  // 1. malloc local variable for the calculation of the sparsity 
  int nz_max_per_row = 150;

  int *col_ldx = (int *)malloc(Msize*nz_max_per_row*sizeof(int));
  int *row_nnz = (int *)malloc(Msize               *sizeof(int));
  for(i = 0; i < Msize; i++) row_nnz[i] = 0;

  int dof_i, dof_j;
  for(elem = 0; elem < N_element; elem++){
    // add each entries to global sparsity pattern
    for(i = 0; i < 6; i++){
      dof_i = ElementEdge[6*elem + i];
      row_start = nz_max_per_row*dof_i;
      for(j = 0; j < 6; j++){
	dof_j = ElementEdge[6*elem + j];
	exist = 0; 
	for(pos = row_start; pos < row_start + row_nnz[dof_i]; pos++){
	  if(col_ldx[pos] == dof_j){
	    exist = 1;
	    break;
	  }
	}
	if(!exist){  // do not exist yet, then add it!
	  col_ldx[pos] = dof_j;
	  row_nnz[dof_i] += 1;
	}
      }
    }
  }

  // 2. compress the sparsity pattern
  *nnz = 0;
  for(i = 0; i < Msize; i++)
    (*nnz) += row_nnz[i];

  (*col_idx) = (int *)malloc((*nnz)*sizeof(int));
  (*row_ptr) = (int *)malloc((Msize+1)*sizeof(int));
  (*row_ptr)[0] = 0;

  // This cycling should be exact!
  // The programmer should responsible for violate memory access!
  for(i = 0; i < Msize; i++){
      for(j = 0; j < row_nnz[i]; j++){
	  // just copy and skip the zero entries
	  (*col_idx)[(*row_ptr)[i] + j] = col_ldx[nz_max_per_row*i + j]; 
	}
      // and set the start pointer of the next row
      (*row_ptr)[i+1] = (*row_ptr)[i] + row_nnz[i]; 
    }

  // 3. free local variable
  free(row_nnz);
  free(col_ldx);
}


/**
 * Sparse pattern with given element dof
 * 
 */
void sparsity_pattern(int **row_ptr, int **col_idx, int Msize, int *nnz,
		      int *ElementDOF, int N_element, int n_dof_per_elem)
{
  int i = 0, j = 0, elem = 0;
  int pos = 0, row_start = 0;
  int exist = 0;

  // 1. malloc local variable for the calculation of the sparsity 
  int nz_max_per_row = 200;

  int *col_ldx = (int *)malloc(Msize*nz_max_per_row*sizeof(int));
  int *row_nnz = (int *)malloc(Msize               *sizeof(int));
  for(i = 0; i < Msize; i++) row_nnz[i] = 0;

  int dof_i, dof_j;
  for(elem = 0; elem < N_element; elem++){
    // add each entries to global sparsity pattern
    for(i = 0; i < n_dof_per_elem; i++){
      dof_i = ElementDOF[n_dof_per_elem*elem + i];
      row_start = nz_max_per_row*dof_i;
      for(j = 0; j < n_dof_per_elem; j++){
	dof_j = ElementDOF[n_dof_per_elem*elem + j];
	exist = 0; 
	for(pos = row_start; pos < row_start + row_nnz[dof_i]; pos++){
	  if(col_ldx[pos] == dof_j){
	    exist = 1;
	    break;
	  }
	}
	if(!exist){  // do not exist yet, then add it!
	  col_ldx[pos] = dof_j;
	  row_nnz[dof_i] += 1;
	}
      }
    }
  }

  // 2. compress the sparsity pattern
  *nnz = 0;
  for(i = 0; i < Msize; i++)
    (*nnz) += row_nnz[i];

  (*col_idx) = (int *)malloc((*nnz)*sizeof(int));
  (*row_ptr) = (int *)malloc((Msize+1)*sizeof(int));
  (*row_ptr)[0] = 0;

  // This cycling should be exact!
  // The programmer should responsible for violate memory access!
  for(i = 0; i < Msize; i++){
      for(j = 0; j < row_nnz[i]; j++){
	  // just copy and skip the zero entries
	  (*col_idx)[(*row_ptr)[i] + j] = col_ldx[nz_max_per_row*i + j]; 
	}
      // and set the start pointer of the next row
      (*row_ptr)[i+1] = (*row_ptr)[i] + row_nnz[i]; 
    }

  // 3. free local variable
  free(row_nnz);
  free(col_ldx);
}

/************************************************************
 only assemb the element-wise linear system into sparse matrix format
*************************************************************/
/**
 *  add entry to sparse matrix
 * 
 */
void add_sparse(int *row_ptr, int *col_idx, 
		double complex *val_nnz, 
		int dof_ii, int dof_jj, 
		double complex entry_val)
{
  int pos;
  for(pos = row_ptr[dof_ii]; pos < row_ptr[dof_ii+1]; pos++){
    if(col_idx[pos] == dof_jj)  break;   // find entry(dof_ii, dof_jj)
  }

  if(pos == row_ptr[dof_ii + 1]){ // wrong in pattern
    printf("Mistake in the pattern of sparse matix.\n");
  }

  val_nnz[pos] += entry_val;    // add the entry value to position       

}

void assemb_mat(int *row_ptr, int *col_idx, double complex *val_nnz,  // CSR formatted
		double complex *element_matrix,          // element-wise sparse matrix
		int *ElementDOF, int N_element, int N_dof_per_elem) 
{
  int e, ii, jj, idx;
  int *DOF;
  double complex *elem_mat;

  struct timeval tv;
  struct timezone tz;
  double t_start, t_end;

  gettimeofday(&tv, &tz);
  t_start = tv.tv_sec + tv.tv_usec*1.0e-06;

  for(e = 0; e < N_element; e++){
    elem_mat = element_matrix + N_dof_per_elem*N_dof_per_elem*e;
    DOF = ElementDOF + N_dof_per_elem*e;

    for(idx = 0, ii = 0; ii < N_dof_per_elem; ii++){
      // parallable according to DOF[ii]

      for(jj = 0; jj < N_dof_per_elem; jj++, idx++){
	add_sparse(row_ptr, col_idx, val_nnz, 
		  DOF[ii], DOF[jj], elem_mat[idx]);
      }
      // end of parallable area

    } // finish element e

  }// run over all elements 

  gettimeofday(&tv, &tz);
  t_end = tv.tv_sec + tv.tv_usec*1.0e-06;
  printf("Spent %f seconds to  assemble the CSR format matrix.\n", 
	 (t_end - t_start));
    
} // end of assemble_APhi_matrix



// only assemb the element-wise rhs into one single vector
void assemb_rhs(double complex *val_rhs, int N_dof,
		double complex *element_rhs,
		int *ElementDOF, int N_element, int N_dof_per_elem)
{
  int e,j;
  /* int *DOF; */
  /* double complex *elem_rhs; */

  for(j = 0; j < N_dof; j++)  
    val_rhs[j] = 0.0 + csqrt(-1)*0.0;

  for(j = 0; j < N_element*N_dof_per_elem; j++)
    val_rhs[ElementDOF[j]] += element_rhs[j];

  /* for(e = 0; e < N_element; e++){ */
  /*   DOF = ElementDOF + e*N_dof_per_elem; */
  /*   elem_rhs = element_rhs + e*N_dof_per_elem; */

  /*   for(j = 0; j < N_dof_per_elem; j++) */
  /*     val_rhs[DOF[j]] += elem_rhs[j]; */
  /* } */

}


/**
 * Apply the Dirichlet Boundary condition
 * 
 */

void apply_Dirichlet(int *row_ptr, int *col_idx,
		     double complex *val_nnz,
		     double complex *val_rhs,
		     int *bnd_flag,
		     double complex *bnd_val,
		     int N_dof)
{
  int i,j,pos,idx;  
  for(i = 0; i < N_dof; i++){

    // 1. skip the non dirichlet nodes
    if(bnd_flag[i] != 1) continue; 

    // 2. set the boundary value
    val_rhs[i] = bnd_val[i]; 

    // 3. Modify the CSR matrix
    for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++){
      j = col_idx[pos];
      if(j == i) {
	val_nnz[pos] = 1.0 + 0.0*csqrt(-1); // diagnal
      }
      else{ // if(cabs(val_nnz[pos]) > 1e-14)
	val_nnz[pos] = 0.0 + 0.0*csqrt(-1); // non-diagnal
	// and column check (j,i), 
	for(idx = row_ptr[j]; idx < row_ptr[j+1]; idx++){
	  if(col_idx[idx] == i){ 
	    // find nnz entry (j,i)
	    val_rhs[j] -= val_nnz[idx]*val_rhs[i]; // reduced to right hand side
	    val_nnz[idx] = 0.0 + 0.0*csqrt(-1); // set to ZERO
	  }
	} // finish clean row j
      } // finish checking (i,j)

    } // finish modify the matrix

  } // finish all dirichlet boundary dof

}


/**
 * Linear solver 
 * calling cholmod sparse decomposition 
 */


void cholmod_solver(int *row_ptr, int *col_idx, double complex *val_nnz, 
		    double complex *val_rhs, double complex *val_sol, 
		    int N_dof)
{
  
}

/**
 * Linear Solver
 * by calling super lu function
 * direct method
 */

void superLU_solver(int *row_ptr, int *col_idx, double complex *val_nnz, 
		    double complex *val_rhs, double complex *val_sol, 
		    int N_dof)
{
  /* SuperMatrix A; */
  /* SuperMatrix L;    */
  /* SuperMatrix U;    */
  /* int      *perm_c; */
  /* int      *perm_r; */
 
  /* int      nrhs, ldx, info; */
  /* SuperMatrix B; */

  /* superlu_options_t options; */
  /* SuperLUStat_t stat; */
  
  /* set_default_options(&options); */
  /* StatInit(&stat); */
    
  /* // 1. convert the complex linear system */
  /* zCreate_CompRow_Matrix(&A, N_dof, N_dof, row_ptr[N_dof],  */
  /* 			 val_nnz, col_idx, row_ptr,  */
  /* 			 SLU_NR, SLU_Z, SLU_GE); */
  /* // rhs */
  /* nrhs   = 1; */
  /* zCreate_Dense_Matrix(&B, N_dof, nrhs, val_rhs, N_dof,  */
  /* 		       SLU_DN, SLU_Z, SLU_GE); */

  /* // 2.solve  */
  /* perm_c = (int *)malloc(N_dof*sizeof(int)); */
  /* perm_r = (int *)malloc(N_dof*sizeof(int)); */
  /* zgssv(&options, &A, perm_c, perm_r, &L, &U, &B, &stat, &info);  */

  /* // 3. set back the value */
  /* int ii; */
  /* //  for(ii = 0; ii < N_dof; ii++) val_sol[ii] = B[ii]; */

  /* // 4. free spaces */
  /* StatFree(&stat); */

  /* SUPERLU_FREE (perm_r); */
  /* SUPERLU_FREE (perm_c); */
  /* Destroy_CompRow_Matrix(&A); */
  /* Destroy_SuperMatrix_Store(&B); */
  /* Destroy_SuperNode_Matrix(&L); */
  /* Destroy_CompRow_Matrix(&U); */
  
}

       
/**
 * print out the sparse matrix when it is small.
 * 
 */

void output_mat(char *fileName, 
		int *row_ptr, int *col_idx, double complex *val_nnz, 
		double complex *val_rhs, 
		int Msize)
{
  int i,j,pos;
  FILE *fout = fopen(fileName, "w");
  int nnz = row_ptr[Msize];

  // 1. size information
  fprintf(fout, "%d \t%d\n", Msize, nnz);

  // 2. sparse matrix
  for(i = 0; i < Msize; i++){  // write line i
    fprintf(fout, "%d\t%d\n", i, (row_ptr[i+1] - row_ptr[i])); // nnz of this row
    for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++){
      fprintf(fout, "\t%d\t %18.8f %18.8f\n", 
	      col_idx[pos],creal(val_nnz[pos]), cimag(val_nnz[pos]));
    }
    fprintf(fout, "\n"); // empty line
  }

  // 3. the rhs
  for(i = 0; i < Msize; i++)
    fprintf(fout, "%18.8f %18.8f\n",creal(val_rhs[i]),cimag(val_rhs[i]));

  fclose(fout);

}


void write_cplin_sys(char *file_name,
		     int *row_ptr, int *col_idx, 
		     double complex *val_nnz, 
		     double complex *rhs,
		     int n_dof)
{
  int i,pos,nnz;
  FILE *fout = fopen(file_name,"w");

  nnz  = row_ptr[n_dof];
  fprintf(fout, "%d %d\n", n_dof, nnz);
  /* for(i = 0; i < N_node; i++) */
  /*   fprintf(fout, "%d %lf\n", row_ptr[i], vec_rhs[i]); */
  /* fprintf(fout, "%d\n", row_ptr[N_node]); */   // the last index or row_ptr
  fprintf(fout, "\n"); //empty line
  for(i = 0; i < n_dof; i++){
    for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++){
      fprintf(fout, "%d %d %lf %lf\n",i, col_idx[pos],
	       creal(val_nnz[pos]), cimag(val_nnz[pos]));
    }
  }
  fprintf(fout, "\n");
  for(i = 0; i < n_dof; i++)
    fprintf(fout, "%f %f\n", creal(rhs[i]), cimag(rhs[i]));
  fprintf(fout, "\n");

  fclose(fout);
}


//////////////////////////////////////////////////
/**
 * add entry to sparse entry
 * 
 */
void add_sparse_int(int *row_ptr, int *col_idx, int *val_nnz, 
		    int dof_ii, int dof_jj, int entry_val){
  int pos;
  for(pos = row_ptr[dof_ii]; pos < row_ptr[dof_ii+1]; pos++){
    if(col_idx[pos] == dof_jj)  break;   // find entry(dof_ii, dof_jj)
  }

  if(pos == row_ptr[dof_ii + 1]){ // wrong in pattern
    printf("Mistake in the pattern of sparse matix.\n");
  }

  val_nnz[pos] += entry_val;    // add the entry value to position       

}  // end


/**
 *  add entry to sparse matrix
 * 
 */
void add_sparse_real(int *row_ptr, int *col_idx, double *val_nnz, 
		     int dof_ii, int dof_jj, double entry_val)
{
  int pos;
  for(pos = row_ptr[dof_ii]; pos < row_ptr[dof_ii+1]; pos++){
    if(col_idx[pos] == dof_jj)  break;   // find entry(dof_ii, dof_jj)
  }

  if(pos == row_ptr[dof_ii + 1]){ // wrong in pattern
    printf("Mistake in the pattern of sparse matix.\n");
  }

  val_nnz[pos] += entry_val;    // add the entry value to position       

}


/**
 *  add entry to sparse matrix
 * 
 */
void add_sparse_complex(int *row_ptr, int *col_idx, 
		     double complex *val_nnz, 
		     int dof_ii, int dof_jj, 
		     double complex entry_val)
{
  int pos;
  for(pos = row_ptr[dof_ii]; pos < row_ptr[dof_ii+1]; pos++){
    if(col_idx[pos] == dof_jj)  break;   // find entry(dof_ii, dof_jj)
  }

  if(pos == row_ptr[dof_ii + 1]){ // wrong in pattern
    printf("Mistake in the pattern of sparse matix.\n");
  }

  val_nnz[pos] += entry_val;    // add the entry value to position       

}


/**
 * get entry value
 * 
 */

double sparse_entry_real(int *row_ptr, int *col_idx, double *val_nnz, int i, int j)
{
  int pos;
  for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++)
    if(col_idx[pos] == j) break;
  if(pos == row_ptr[i+1])
    return 0.0; // (i,j) are not int the pattern, so just return 0.0
  else
    return val_nnz[pos];  
}


/**
 * get entry value
 * 
 */

double complex sparse_entry_complex(int *row_ptr, int *col_idx, 
				    double complex *val_nnz, 
				    int i, int j)
{
  int pos;
  for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++)
    if(col_idx[pos] == j) break;
  if(pos == row_ptr[i+1])
    return 0.0; // (i,j) are not int the pattern, so just return 0.0
  else
    return val_nnz[pos];  
}

/**
 * Apply the dirichlet boundary to the sparse linear system.
 * Please note that this is pure linear operation,
 * user is responsible for calculate the specific value in their main routine
 */

void apply_const_complex(int dim, int *row_ptr, int *col_idx, 
			 double complex *val_nnz, double complex *rhs, 
			 int nn, int *dof_indices,double complex *val)
{
  int ix, ik, ii, jj;
  int dof_ii, dof_jj;
  
  for(ix = 0; ix < nn; ix ++){

    dof_ii      = dof_indices[ix]; 
    rhs[dof_ii] = val[ix];  // set the value
        
    for(ik = row_ptr[dof_ii]; ik < row_ptr[dof_ii + 1]; ik++){ 
      dof_jj = col_idx[ik];
      if(dof_jj == dof_ii){    // elliminate line dof_ii
	val_nnz[ik] = 1; 
      }
      else{
	val_nnz[ik] = 0; 
	for(jj = row_ptr[dof_jj]; jj < row_ptr[dof_jj+1]; jj++){
	  if(col_idx[jj] == dof_ii){  // find (dof_jj, dof_ii)
	    rhs[dof_jj] -= val_nnz[jj]*rhs[dof_ii];  // put it to rhs
	    val_nnz[jj] = 0;   // and elliminate it in the matrix
	    break;   
	  }
	} // finish (dof_jj,dof_ii) into rhs  
      }  // clear non-diagnal entries
    }

  }  // finish treating all nnz in the row dof_ii 
 
}


/*
  used to check the symmetricity of the sparse matrix
 */
int test_symm(int *row_ptr, int *col_idx, int *val_nnz, int Msize)
{
  int i, j, pos, pos_symm;
  int symm = 1;
  double dist = 0;
  for(i = 0; i < Msize; i++){
    for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++){
      // check each nnz

      j = col_idx[pos];
      for(pos_symm = row_ptr[j]; pos_symm < row_ptr[j+1]; pos_symm++){
	if(i == col_idx[pos_symm]) break;   // find the symmetric element
      }

      if(pos_symm == row_ptr[j+1]) 
	printf("Even the matrix pattern is not symmetric!\n");

      dist = abs(val_nnz[pos] - val_nnz[pos_symm])*1.0;
      if(dist > 1e-3){ // check symmetric via relative error
	symm = 0;
	break;
      }

    } // finish row i
  }  // finish all nnz

  return symm;
}

/**
 * write out the sparse matrix nnz into txt files
 * 
 */

void write_txt(int *row_ptr, int *col_idx, double complex *val_nnz, int Msize, char* fileName)
{
  FILE *fout = fopen(fileName, "w");
  fprintf(fout, "Msize nnz\n");
  int i,pos;
  for(i = 0; i < Msize; i++){  // write line i
    fprintf(fout, "%d\t%d\n", i,(row_ptr[i+1] - row_ptr[i])); // nnz of this row
    for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++){
      fprintf(fout, "\t%d\t %18.8f + %18.8f i\n", 
	      col_idx[pos],creal(val_nnz[pos]), cimag(val_nnz[pos]));
    }
    fprintf(fout, "\n"); // empty line
  }

  fclose(fout);

}

// the second version with different data structure
void write_txt2(int *col_idx, double complex *val_nnz, int *row_size,
		int Msize, int nz_row, char* fileName)
{
  FILE *fout = fopen(fileName, "w");
  fprintf(fout, "Msize nnz\n");
  int i,j,p_start;
  for(i = 0; i < Msize; i++){  // write line i
    fprintf(fout, "%d\t%d\n", i, row_size[i]); // nnz of this row
    p_start = i*nz_row;
    for(j = 0; j < row_size[i]; j++){
      fprintf(fout, "\t%d\t %18.8f + %18.8f i\n", col_idx[p_start+j],
	      creal(val_nnz[p_start+j]), cimag(val_nnz[p_start+j]));
    }
    fprintf(fout, "\n"); // empty line
  }

  fclose(fout);

}
